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ABSTRACT 

We use the same physical model to simulate four galaxies that match the relation 
between stellar and total mass, over a mass range that includes the vast majority of 
disc galaxies. The resultant galaxies, part of the Making Galaxies in a Cosmologi- 
cal Context (MaGICC) program, also match observed relations between luminosity, 
rotation velocity, size, colour, star formation rate, HI mass, baryonic mass, and metal- 
licity. Radiation energy feedback from massive stars and supernova energy balance the 
complex interplay between cooling gas, regulated star formation, large scale outflows, 
and recycling of gas in a manner which correctly scales with the mass of the galaxy. 
Outflows, driven by the expansion of shells and superbubbles of overlapping super- 
nova explosions, also play a key role in simulating galaxies with exponential surface 
brightness profiles, flat rotation curves and dark matter cores. Our study implies that 
large-scale outflows are the primary driver of the dependence of disc galaxy properties 
on mass. We show that the degree of outflows invoked in our model is required to meet 
the constraints provided by observations of OVI absorption lines in the circum-galactic 
media of nearby galaxies. 
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1 INTRODUCTION 



The rotation velocity, size, lumi nosity (e.g. Courteau et al.l 
l2007h , st ar formation rat e (e.g. ISalim et al.l l2007h . stellar 
mass fe .g.lBell et al.1 2003). total viri al mass (incl uding dark 
matter, Mandelbau m et al. | e.g. 120061 : iMore et al.1 e.g. l201lh. 
of neutral hydrogen gas (e.g. IVerheiien &; Sancisil 



|200l|), total baryon mass fe.g. lMcGaughll2005h and the abun- 

Till 



danc e of oxygen relative to hydrogen fe.g. lTremonti et al.l 
120041 ) are observationally well-determined. Together, the re- 
lation between these properties provide stringent constraints 
which simulations of galaxy formation must satisfy. 

A large number of such disc galaxy properties vary 
with mass. Growing evidence suggests that the ratio of stel- 
lar mass to total mass of galaxies is low compared to the 
ratio of baryons to dark matter in the Universe. The re- 
lation is steep in the mass range where most disc galax- 
ies exist, below a stellar mass ~2xlO lo M0, where a dou- 
bling of halo mass results in a factor of eight more stellar 
mass. The relation has been determined by rank-ordering 
the masses of dark matter halos formed in simulations within 
the cold dark matter cosmological framework, and then 



rank-ordering galaxies based on their stellar mass as deter- 
mined in large scale galaxy surveys, and directly matching 
stellar and halo masses , assuming a monotonic relati on be- 
tween the two ([Moster et al.1 l2010l : iGuo et al.1 boioh . The 
results agree remarkably well with direct derivations of the 
relation, where dark matter halo masses are d etermined by 
gravitational lensing ([Mand elbaum et al. 2006) and the dy- 
namics of satellite galaxies ([More et al.l 1201 lh . 

Recent progress in simulating galaxies in a cosmological 
cold dark matter framework has resulted in the ability to 
match a n increasing number of properties of observed disk 
galaxies (iGoyernato et al.1 2010l: IPiontek &; Steinmet j|201ll : 



lAgertz et alj|201ll : iGuedes et al.l l201lh but have not been 
able to model the dependence of galaxy properties on mass. 

Two approaches have emerged to regulate star forma- 
tion in simulations. One meth od is to lower star forma- 
tion efficiency either directly (lAgertz et al.l 201 ID or us- 
ing small scale star formation physics based on the avail- 
ability of molecular hydrogen (H2 ) from which stars form 
([Kuhlen et al.ll201lh . lAgertz et all ([201 lh produced a sim- 
ulated galaxy which shares ma ny features with the Milky 
Way, while ([Kuhlen et al.ll201lh produced low mass galaxies 
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at high redshift that match empirical star formation laws. 
Yet these models provide no mechanism to account for the 
steep decrease in stellar mass of galaxies found in lower mass 
galaxies, even accounting for the fact that lower mass galax - 
ies have higher gas contents (e.g. IPeeples &; ShankarluOllh . 

A contrasting model invoking energy "feedback" from 
supernovae to regulate star formation has successfully sim- 
ulated dwarf and Milky Way mas s galaxies that share 
many properties of ob s erved galaxies ( Govern ato et al.ll2010l : 
IPiontek fe Steinmetzl 1201 ll : fauedes et al.1 l201lh . Despite 
these successes, the galaxies simulated so far with the super- 
nova energy feedb ack models form too many stars relative 
to their total mass (ISawala et al.|[201ll : IPiontek &; Steinmetzl 
[2Qlll: lAvila-Reese et al.l l201lTl7 particularly in low mass sim- 
ulations. Additionally, the simulations fail to matc h the spe- 
cific star formatio n rates of observed galaxies (|Colm et al.1 
120101 : lAvila-Reese et al.ll201lh . 

In this paper we take the latter approach of directly 
injecting SNe energy. A clue to improving on such mod- 
els came from measuring the velocity dispersion of the gas 
of our simulated dwarf galaxies, which was found to be 
significantly less turbu lent than observed dwarf galaxies 
(Pil kington et al.ll201lh . suggesting that "feedback" energy 
has been under-estimated. Recent evidence suggests that ra- 
diative feedback from massive stars before they explode as 
supernovae can have significant effects in regulating star for- 
mation and adding e nergy to gas that surrounds newly form- 
ing clusters of stars (iNath fe Silkll2009l : iMurrav et al] 1201 it 
iHopkins et al.ll201lLl2012h . Young, massive stars exert forces 
from their UV photons, stellar winds, and warm gas pressure 
from photo-ionized regions. In massive galaxies, momentum 
driven winds due to radiation pressure may be the primar y 
driver of outflows ([Murray et al.ll2005l : [Hopkins et al.ll2012l ). 
However, in galaxies of mass relevant to our study, the effect 
of the radiation feedback from massive stars is primarily to 
(i) prevent collapse of gas to small, dense regions, dissociat- 
ing GMCs and regulating star formation, and (ii) stir up gas 
in star forming regions, punching hole s that allow SNe rem - 
nants to expand and drive outflows ([Hopkins et al.l [2Q12h , 
and (iii) provide pressure support in the disk. 

Our model does not resolve this local radiation feed- 
back from star forming regions nor individual supernovae, 
but it does resolve the shells and bubbles that are created 
by overlapping supernovae. Without the ability to resolve 
the details, we employ a relatively crude thermal implemen- 
tation of radiation feedback from massive stars, with the 
aim being to mimic their most important effects on scales 
that we resolve, i.e. to regulate star formation, and enhance 
inhomogeneity in the ISM, creating turbulence that better 
matches observations, and to allow the expansion of the SNe 
driven super-bubbles which drive outflows. 

Using this model, we simulate four galaxies with stel- 
lar masses ranging from 2xl0 8 to 1.4xlO lo M0. Crucially, 
the simulated galaxies (i) were run at the same resolution 
with identical input physics, and (ii) vary in their mass as- 
sembly history. Comparing the two most massive galaxies, 
SG3 has most of its merger activity at early times, prior to 
z = 2, whilst in SG4 significant merging activity occurs after 
z = 1. To illustrate the effect of feedback, we include results 
of SGI using less feedback (SG1LF). It forms a galaxy with 
too many stars and is too compact. To study the way our 
model varies with resolution, we re-simulated SG3 at 8 times 



lower resolution (SG1LR), run with the same feedback pre- 
scriptions, as well as a more massive galaxy SG5LR which 
extends the mass range of our sample out to 3xlO lo M0. 
These lower resolution simulated galaxies also match the 
scaling relations. 

Section [2] describes the code and the initial conditions. 
Section [3] presents the basic properties of the simulations in- 
cluding star formation histories, rotation curves, radial light 
and dark matter density profiles. In Section U we show that 
the galaxies all fit on the observed relations between rotation 
velocity, size, luminosity, star formation rate, stellar mass, 
virial mass, HI mass, total baryon mass and metallicity. Sec- 
tion [5] shows that the outflows, which are central to our mod- 
els, produce metal-enriched hot halos that reproduce the ob- 
served OVI column densities in the circum-galact i c-me dium 
(CGM, IProchaska et al.l l201ll ; iTumlinson et all I Hoill)- In 
Section [6] we summarise our model for disc galaxy forma- 
tion. 



2 THE SIMULATIONS 

2.1 The simulation code: GASOLINE 



We have used GASOLINE ([Wadslev et al.ll2004h . a fully paral- 
lel, N-Body+smoothed particle hydrodynamics (SPH) code, 
to compute the evolution of the collisionless and dissipative 
components in the simulations. The essential features of the 
code are outlined here. 

Fluid elements representing gas are inte- 
grated using Smooth Particle Hydrodyn amics (SPH, 
iGingold fe Monaghan! Il977l : iMonaghanl Il992h . GASOLINE is 
fully Lagrangian, spatially and temporally adaptive, and 
efficient for large N. Dissipation in shocks is modelled using 
the quad ratic term of the standard Monaghan artificial 
viscosity (Monaghan 1992). The Balsara correction term 
is used to reduce unwanted shear viscosity. GASOLINE 
employs cooling d ue to H, He, and a variety of metal lines 
( She n et "all l2010h . The metal cooling grid is con structed 
using CLOUDY (version 07.02 iFerland et alJ[l998 ). assum- 
ing ionisation equilibrium. A unifor m ultraviolet ionising 
background ([Haardt &; Madaul 1 19961 ) is used in order to 
calculate the metal cooling rates self-consistently. Two 
mechanisms are used to prevent gas from collapsing to 
higher densities than SPH can physically resolve: (i) to 
ensure that gas resolves the Jeans mass and does not 
artificially fragment, pressure is added to the gas i n high 
density star forming regions ([Robertson &; Kravtsovll2008h , 
(ii) a maximum density limit is imposed by setting a 
minimum SPH smoothing length of 0.25 times that of the 
gravitational softening length of e = 155pc. Each simulation 
has between 5-15 million particles within the virial radius 
at z = 0, with mean stellar particle mass of ~4800 Mq 
( ~ 38000 M© in the low resolution runs). 



2.2 The Initial Conditions 

The simulations described here are cosmological zoom simu- 
lations derived f rom the McMaster Unbiased Galaxy Simu- 
lations (MUGS) (IStinson et alJ lioioh . Here, those initial con- 
ditions are scaled down, so that rather than residing in a 68 
Mpc cube, it is inside a cube with 34 Mpc sides. This resizing 
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Table 1. Simulation data 



Name MUGS gas part. IMF c* M halo M* M R S* n** 

label mass [M ] [M ] [M ] 



SGI 


g5664 


2.5xl0 4 


Chabrier 


0.17 


6.5xl0 10 


2.3xl0 8 


-17.0 


21.2 


1.2 


0.82 


SG2 


gl536 


2.5xl0 4 


Chabrier 


0.17 


8.3xl0 10 


4.5xl0 8 


-17.5 


21.4 


1.4 


0.77 


SG3 


gl5784 


2.5xl0 4 


Chabrier 


0.17 


1.8xlO n 


4.2xl0 9 


-20.0 


20.0 


2.0 


1.17 


SG4 


gl57807 


2.5xl0 4 


Chabrier 


0.17 


3.2xlO n 


1.4xl0 10 


-21.2 


19.0 


2.5 


1.23 


SG1LF 


g5664 


2.5xl0 4 


Kroupa 


0.17 


7.3xl0 10 


8.7xl0 9 


-19.9 


18.9 


1.0 


1.4 


SG3LR 


gl5784 


2.0xl0 5 


Chabrier 


0.33 


1.8xlO n 


4.6xl0 9 


-20.1 


19.9 


2.1 


1.1 


SG5LR 


g!536 


2.0xl0 5 


Chabrier 


0.33 


7.6xlO n 


3.0xl0 10 


-21.7 


20.5 


3.6 


1.3 



* central surface brightness (/io), scale-lengths (size, S) and sersic indices (n) from fits in the I band 
**SGl-4 and SGLR are single sersic fits, SG1LF and SG5LR are 2 component sersic bulge+exponential disk fits 



allows us to compare galaxies with exactly the same merger 
histories at a variety of masses. Differences in the underly- 
ing power spectrum that result from this re s caling are mi- 
nor (Spri ngel et al . 2008; Maccio et al.l [20081 ; iKannan et al.1 
I2012h 7 and do not effect our results. SG5LR is run with the 
original 68 Mpc cube, having the same initial conditions as 
SG3. 

Table 1 shows the properties of the galaxies that were 
selected and rescaled from the MUGS sampl e. Note that 
SG3 uses the same initial conditions used in Brook et al.l 
<l2012h . which showed the impact of galactic fountains on an- 
gular momentum distribution of galaxies. Howe ver, SG3 uses 
a dif ferent initia l mass function (IM F) from ([Brook et al.l 
l2012h . repl acing [Kroupa et all ([1993) with the more com- 
monly used Chab rierl d2003l) . which is more top heavy. The 
conclusions of Brook et al.l ([20121 ) remain valid in the new 
simulations. 



2.3 Star Formation and Feedback 

When gas reaches cool temperatures (T < 10, 000 K) in a 
dense environment (n t h > 9.3 cm -3 ), it becomes eligible to 
form stars. This value for nth is the maximum density gas 
can reach using gravity, 32 m^ as /e 3 . Such gas is converted 
to stars according to the Schmidt Law: 



AM, 
At 



tdy 



(1) 



where AM* is the mass of the stars formed in At, the time 
between star formation events (0.8 Myrs in these simula- 
tions), m gas is the mass of the gas particle, tdyn is the gas 
particle's dynamical time and c* is the efficiency of star for- 
mation, in other words, the fraction of gas that will be con- 
verted into stars during tdyn- Effective star formation rates 
are determined by the combination and interplay of c* and 
feedback, and so degeneracies do exist between feedback en- 
ergy and the value of c*. In this study, c* is ultimately the 
free parameter that sets the balance of the baryon cycle off 
cooling gas, star formation, and gas heating. In our fiducial 
runs, c*=1.67%. 

Stars feed energy back into surrounding gas. Two types 
of energetic feedback are considered in these simulations, su- 
pernovae and early stellar radiation feedback from massive 
stars. Supernova feedback is implem ented using the blast- 
wave formalism ([Stinson et al.l l2006) and deposits 10 51 erg 
of energy into the surrounding medium at the end of the 
stellar lifetime of every star more massive than 8 M©. Since 



stars form from dense gas, the energy would be quickly ra- 
diated away due to the efficient cooling. For this reason, 
cooling is disabled for particles inside a blast region of size 
R = lO 1 - 74 ^! 32 ^ 016 P 4°- 20 pc and for the length of time 
t = 10 6 - 85 ^5 i 32 n8- 34 Pn" 4 - 70 yr. Here, E 51 = 10 51 



erg, n is 



the ambient hydrogen density, and P04 = 10 _ Pok~ , where 
Po is the ambient pressure and k is the Boltzmann constant. 
Both no and Po are calculated using the SPH kernel for the 
gas particles surrounding the star. 

Metals are ejected from type II supernovae (SNII), type 
la supernovae (SNIa), and the stellar winds driven from 
asymptotic giant branch (AGB) stars. Ejected mass and 
metals are distributed to the nea rest neighbour gas p arti- 
cles using the smoothing kernel Jst inson et al.l l2006f) , us- 
ing lit erature yields for SN II ([Wooslev &; Weaver! 1995) and 
SNIa ([Nomoto et al.lll99<t ). Metal diffusion is also included, 
such that unresolved turbul ent mixing is trea ted as a shear- 
dependent diffusion term ( Sh en et "all [2010) . This allows 
proximate gas particles to mix their metals. Metal cooling 
is calculated based on the diffused metals. 

Radiation energy feedback from massive stars has been 
included in our model. To model the luminosity of stars, 
a simple fit of the mas s- lumin o sity r elationship observed in 
binary star systems bv lTorred (2010) is used: 



L 



( M \4 



M < 10M Q 
M > 10M Q 



(2) 



Typically, this relationship leads to 2 x 10 50 ergs of energy 
being released from the high mass stars per M© of the entire 
stellar population over the 4.5 Myr between the star's for- 
mation and the commencing of SNII. These ph otons do not 
coup le efficiently with the surrounding ISM ([Frever et al.l 
2006). We thus do not want to couple all of this energy to 
the surrounding gas in the simulation. To mimic this highly 
inefficient energy coupling, we inject 10% of the energy as 
thermal energy in the surrounding gas, and cooling is not 
turned off for this form of energy input. It is well established 
that such thermal energy injection is highly inefficient at the 
spatial and temporal r esolution of the type of cosm ological 
simulations used here (lKatzlll992l : iKav et alJl2002h . This is 
primarily due to the characteristic cooling timescales in the 
star forming regions being lower than the dynamical time. 

In star forming regions where gas has density of ~ 10 
cm -3 , early stellar feedback typically heats the surround- 
ing gas to between ~ 1 to a few xl0 6 K. At the density 
where stars form in our simulations, i.e. n~ 10 cm -3 , cool- 
ing times are significantly shorter than dynamical times, 
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t 1 

50 kpc 

Figure 1. Mock images of the four simulated galaxies, made with SUNRISE dJoussonll2006h that uses the ages and metallicities of simulated 
star particles, and the effects of dust reprocessing, to create stellar energy distributions. Filters which mimic Sloan Digital Sky Survey 
(SDSS) bands i, r, and g, are assigned to colours R, G, and B respectively to produce these images. The galaxies, from left to right, SGI, 
SG2, SG3 and SG4, are shown face-on (top panels) and edge-on (bottom panels). The significant asymmetry in SG4 is caused by the 
late accretion of a low mass satellite galaxy. 



even for temperatures of o rder 10 6 K (see e.g. Figure 1 in 
iDalla Vecchia &: Schavell2012h , typical of those which to our 
gas is heated by our radiation energy feedback, meaning 
that it is only gas particles which escape to less dense re- 
gions that have any affect from this feedback. We found that 
only ~10% of heated particles had dynamical times longer 
than cooling times at high redshift (3>z>l) when the ISM 
is most turbulent and the metallicity is low, reducing to only 
~3% by 2; = 0. Thus, between 90 — 97% is radiated away 
within a single dynamical time, meaning that our effective 
efficiency of coupling radiation energy feedback to the ISM 
is between 0.3—1%. It also means that our implementation 
does not evenly heat gas, i.e. we do not effectively couple 
the <1% of energy to all gas particles affected by radiation 
energy. Rather, our scheme is essentially stochastic, with a 
small number of gas particles effected. This scheme is thus 
somewhat similar to the SNe feedback implementation of 
Dalla Vecch ia fc Schavd (2012). This radiation energy feed- 
back does not directly drive outflows: implemented in the 
absence of SNe feedback results in no outflows. Rather, it 
helps to regulate star formation, and enhances inhomogene- 
ity in the ISM, creating increased turbulence, and allowing 
the expansion of the SNe driven super-bubbles which drive 
outflows. 

To understand the effect of the strengthened feedback, 
a version of SGI was simulated using the stellar feedback 
as implemented in IStinson et all fcOlOh so that stars form 
when gas reached 1 cm -3 , with a lKroupa et all (|l993h IMF, 
with 0.4 xlO 51 ergs deposited per supernova explosion. We 



refer to this as the lower feedback model (LF) in this study, 
but note that this feedback strength is comparable to most 
implementations that a re currently run in the literature 
([Scannapieco et aT]|201lh . 



2.4 Resolution Dependence 

To study the resolution dependence of our model, we also 
simulated two galaxies at 8x lower resolution: SG3LR and 
SG5LR. Clearly, SG3LR is the low resolution version of 
SG3, while SG5LR is a galaxy with stellar (total) mass 
of 3.0xl0 10 (7.6xl0 n )M o . As with other galaxy forma- 
tion simulations in the literature, galaxy properties are not 
precisely the same at diffe rent resolutions when t he same 
parameters are used (e.g. IScannapieco et aD l201lh . In this 
study we aim to retain the same baryon cycle at the two dif- 
ferent resolutions as this drives the simulated galaxy proper- 
ties. Hence we want our star formation rates to match at the 
two resolutions. To achieve this we adjust our free param- 
eter c*, the input star formation efficiency, to ensure that 
the star formation rate remains the same at the two resolu- 
tions. With c*=0.033, SG3LR has a star formation rate that 
closely matches our fiducial high resolution run, SG3. 

As the feedback implementation is the same between 
the two resolutions, the resultant galaxies have the same 
balance of star formation and feedback, and hence very sim- 
ilar baryon cycles, meaning that same key processes occur 
on the scales that we resolve i.e. the same the baryon cy- 
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Figure 2. The star formation histories of the four simulated 
galaxies, SGI (cyan line), SG2 (yellow), SG3 (blue), and SG4 
(red). SGI and SG2 have a bursty, yet reasonably consistent star 
formation histories. SG3 has a peak at high redshift, when it has 
its merging epoch, while SG4 has a much later merging epoch, 
resulting in a star formation rate that peaks at a later time. 



cle between star formation, hot and cold gas, outflows and 
the density of gas in the star forming regions. At z=0, total 
baryon mass within the virial radius, the mass of stars, mass 
of hot gas (>4xl0 4 K), and cold gas (<4xl0 4 K) in SG3LR 
are all within 5-10% of the values for SG3. We will see that 
the resultant galaxies have almost identical properties with 
respect to scaling relations. We subsequently use the same 
calibration and input parameters for SG5LR, a more mas- 
sive galaxy. This shows that we can form reasonable galaxies 
even at relatively low resolution, and it allows us to increase 
the range of masses over which our simulated galaxies match 
the scaling relations. 
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Figure 3. The growth of the virial mass of the halos of the sim- 
ulated galaxies. The later build up of mass in SG4, which is still 
undergoing significant merging activity well after z=l, relates di- 
rectly to its later peak in star formation. 



3.2 Surface Brightness Profiles 

Figure U shows exponential (blue line) and single sersic (red 
line) fits to the I band profiles of the face-on surface bright- 
ness maps of each simulated galaxy, measured from face-on 
images and including the effects of dust reprocessing using 
SUNRISE. In each case we fit from the centre out to the radius 
where the I band surface brightness is 25Marcsec -2 . The 
central surface brightness (/x ), exponential scale- length (S), 
and sersic index (n) are shown. The values of scale-length (S) 
are used for the "size" of the galaxy in this paper, to make 
consistent wi th the observational d ata with which compari- 
son is made ([Courteau et al.l [2QQ7h . Each galaxy exhibits a 
nearly exponential surface brightness profile. In every case, 
the Sersic index, n < 1.5. 



3 SIMULATED GALAXY CHARACTERISTICS 

Mock observations of the four main simulated galaxies 
are shown in Figure [1] These images were creat ed using 
the M onte Carlo radiative transfer code SUNRISE ([jonssonl 
2006), which simulates the effect of dust absorption. The 
morphologies of the galaxies range from dwarf irregular 
(SGI, SG2) to more structured disc galaxies (SG3, SG4). 
The significant asymmetry in SG4 is caused by the late ac- 
cretion of a low mass satellite galaxy. 

3.1 Star Formation and Mass Accretion 

Figure [2] shows the star formation histories of the simulated 
galaxies, SGI (cyan line), SG2 (yellow), SG3 (blue) and SG4 
(red). SGI and SG2 both have bursty, yet reasonably con- 
sistent star formation histories. SG3 has a peak at high red- 
shift, when it has its merging epoch, while SG4 has a much 
later merging epoch, resulting a star formation rate that 
peaks at a later time. The growth of the virial mass of the 
halos of the simulated galaxies is shown in Figure [3] Merger 
events appear as significant increases in the total mass. The 
later build-up of mass in SG4 relates directly to its later 
peak in star formation. 



3.3 Rotation Curves 

Figure [5] shows the rotation curves of four main galax- 
ies along with the low feedback version of SGI. The rota- 
tion velocities (V c ) are computed by calculating the mass 
within spherical radial shells M(r), then relating via V c (r) = 
Y (GM(r)/r) where G is Newton's gravitational constant. 
Radial units are scaled by the disc scale- length in each case. 
The excessive central mass distr ibution that has p lagued 
simulations of galaxy formation (|Maver et al.l [2008h is ab- 
sent in each case, reflected in the lack of central "peak" in 
the rotation curves. The vertical dotted and dashed lines are 
at 2.2 and 3.5 scale- lengths, where observati onal determina- 
tions of ro tational velocity are often made (|Courteaulll997l : 
iGiovanelli et alJll997h . We can see that in SGI, the rota- 
tion curve is still rising quite rapidly at 2.2 scale-lengths. 
This issue of rising rotation cu rves is also commo n in ob- 
servations of low mass galaxies (|Epinat et all [2008). In this 
paper, we report the value of V c at 3.5 scale-lengths in each 
case, as it best reflects the underlying mass of the system. 
Also shown is SG1LF, as a dashed cyan line. The contrast 
between the two SGI galaxies is extreme, with much of the 
available baryons forming stars in the central region of the 
low feedback run, making the central region very dense and 
the rotation velocity very high in this region. 




Figure 4. The I band surface brightness profiles of SGI, SG2, SG3 and SG4 from left to right respectively. Fits are made using both 
exponential (blue lines) and sersic (red lines) profiles. Central surface brightness (/x ) and scale-lengths (S) of the exponential fits as well 
as the sersic index (n) from the sersic fits are noted in the bottom left of each panel. In each case fits are made from radius=0 out to 
I=25Marcsec -2 . 



3.4 Dark Matter Profiles: Cores not Cusps 

The dark matter density profile, on a log-log scale, is plot- 
ted for SGI (cyan line), SG2 (yellow), SG3 (blue), and SG4 
(red) in Figure [6] Also plotted is a dark matter only sim- 
ulation of SG3, which has the steep inner profile that i s 
characteristic of dark matter halos ([Navarro et al.l [l995). 
Radial units are scaled by the virial radius in each case. 
In all the simulated galaxies, the dark matter profile is less 
steep than the profiles in pure dark matter simulations. This 
flattening of the central dark matter profile when baryons 
are included has also been found in cosmolo gical simula- 
tions of dwar f galaxies ([Goyernato et al.ll2010h and massive 
disc galaxies (iMaccio et al.ll2012h . This contrasts with sim- 
ulations that do not include strong feedback, which invari- 
ably have dark matter profiles that are even steeper than 
in pure dark matter simulations ([Gnedin et al lHoill) be- 
cause the dark matter adiabatically contracts as gas cools 
to the center. SG1LF (dashed cyan line) shows exactly this 
behaviour as its dark matter density profile is significantly 
steeper than the simulation that only includes dark matter. 
The fact that the amount of the outflows from the central 
region in the present study results in matching the scal- 
ing relations, favours the findings of simulations where "ex- 
pansion" rather than contraction occurs over the masses of 
these simulations, as semi-anal ytic models have predicted 
(jDutton fe van den Boschll2009h . 



4 SCALING RELATIONS 

The crucial test of galaxy formation models does not 
come from matching a single relation, as that can often 
be achieved through a parameter search of input physics. 
Rather, models are judged by their ability to simultaneously 
match a range of observational constraints. 

4.1 Definitions 

Halo Mass: To make a direct comparison to the abundance 
matching results, w e use M200 from our dark matter only 
runs. As described in lStinson et al ] (|2010l ). one step in creat- 
ing the initial conditions for the MUGS galaxies was running 
a zoom dark matter only simulation at the same resolution 
as the simulation that includes baryons. M200 is the mass 



enclosed inside r2oo , the radius at which the density reaches 
200p c , where p c is the critical density. 

Stellar Mass: The sum of the mass of all stars within the 
radius where the I band surface brightness is 25 M arcsec -2 . 
HI mass: The Saha equation is solved to determine an ion- 
ization equilibrium. This remains an approximation since 
an accurate model of HI mass would require full radiative 
transfer included in the code and is beyond the scope of this 
paper. In particular self shielding from the UV background 
is not included in our model and may result in our derived 
HI masses being under estimates, while photo-ionization of 
HI from the galaxy itself it also excluded. 
Cold gas mass: 4/3 xM///(5/i r ). That is 4/3 the HI mass 
contained within five scale-lengths of the galaxy. A factor 
between 1.2 and 1.4 is commonly used in observations to ac- 
count for heavy elements. We use 4/3 as it is used in the rel- 
evan t observational studies to which we compare ([McGaughl 
l20Q5h . 

Magnitudes, luminosities and colours: come directly from 
the SUNRISE outputs, and include the affects of dust repro- 
cessing on the spectral energy distribution. 



4.2 Relationships 

Panel (a) of Figure [Ushows the relationship between the stel- 
lar mass and host halo mass of the four fiducial simulated 
galaxies as yellow diamonds, over-plotted on the published 
iGuo et al . (2010) relationship. Our free parameter, c*, was 
tuned to place SG3 on this relationship. There was no guar- 
antee that the other simulated galaxies should fall on this 
relationship. The red triangle represents the low feedback 
case, SG1LF, which demonstrates that low feedback results 
in more than an order of magnitude too many stars. 

The green diamonds represents the low resolution runs. 
As stated, we adjusted c* such that the star formation rate 
in SG3LR matched as closely as possible to that of SG3. 
As we have not changed our feedback recipe, the balance 
between star formation and energy feedback remains the 
same at the two resolutions. What is interesting is that using 
these parameters for a more massive galaxy, SG5LR, the 
simulation again fits on the relation, extending the mass 
range over which our balance between star formation and 
feedback results in the correct scaling of stellar mass with 
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Figure 5. Rotation velocity (V c ) versus radius of SGI (cyan 
line), SG2 (yellow), SG3 (blue) and SG4 (red). Radial units are 
scaled by the disc scale-length in each case. The low feedback 
simulation, SG1LF (dashed cyan line) shows the extreme differ- 
ence in stars formed compared to the fiducial run, particularly in 
the central regions, combined with adiabatic contraction rather 
than expansion. The dotted and dashed vertical lines are at 2.2 
and 3.5 scale-lengths. In the paper, we use V c as measured at 
3.5 scale-lengths in each case, as it better reflects the underlying 
mass of the system. 

halo mass. A more thorough analysis of SG5LR will be made 
in a forthcoming study. 

Figures[7l&[8] show the relationships between the rota- 
tion velocity, size, luminosity, star formation rate stellar 
mass, total mass, HI mass, total baryon mass and metallic- 
ity of the simulated disc galaxies. The fiducial simulations 
match the relations remarkably, and their properties cor- 
rectly scale over the mass range of the simulations. Again, we 
emphasise that none of the parameters were tuned in order 
to match any of these relations. As stated, the only tuning 
was to adjust c* such that SG3 (and SG3LR) fits on the stel- 
lar mass-halo mass relation. The low resolution cases (green 
diamonds) also match the relations, with SG3LR matching 
well its high resolution counterpart. 

It is very interesting that the thermal feedback imple- 
mentation results in stellar mass-halo mass relation that has 
such steep dependence on mass (panel 7a). Our feedback has 
no mass dependence nor mass loading included in its imple- 
mentation. It simply inputs a given amount of energy per 
unit of star formation, resulting in pressure from thermal en- 
ergy that drives outflows, controlling the supply of cold gas 
available for star formation. Feedback processes then further 
regulate the star formation rates of this cold gas. 

The other panels in Figure 7 (b-d) relate to the angular 
momentum of the galaxies at given luminosities. Matching 
the Tully-Fisher relation over a wide range in luminosities, 
combined with making disk galaxies of the correct size (with- 
out large stellar bulges) , indicates that we have attained the 
correct amount and distribution of angular momentum. The 
angular momentum of the baryons in these simulations do 
not simply follow the dark matter, as analytic models of- 
ten assume. Rather, the distributions arise from a cycle of 
cooling to the star forming regions, o utflows, and galactic 
fountains fsee llBrook et al.ll201ll . l2012h . 

The relations in Figure 8 put further constraints on the 
baryon cycle of galaxies of varying mass. Even though our 
feedback drives outflows, it does not simply eject all the cold 
gas from the disk. This is the result of the inhomogeneous in- 



Figure 6. The dark matter density profiles, on a log-log scale, 
where radial units are scaled by the virial radius. Also plotted 
is a dark matter only simulation of SG3, which has the steep 
inner profile that is characteristic of dark matter halos. The low 
feedback simulation SG1LF (dashed cyan line) shows significant 
adiabatic contraction. 

terstellar medium, with pockets of hot gas building pressure 
and driving outflows. This is key in allowing the high HI frac- 
tions in the low mass galaxies (Panel 8c) , even though they 
lose a high fraction of their baryons in the outflows. These 
processes also relate directly to the success of the low mass 
simulations in matching the specific star formation rates. 
Simulations have largely had trouble mat ching the specific 
star formation rate for l ow mass galaxies (|Colm et al.ll201ol ; 
lAvila- Reese et all l2Qllh . With low feedback, the star for- 
mation rates follow the gas accretion rates, which exponen- 
tially decline after z ~ 2, meaning that the stars form too 
early. Attempts to form significantly fewer stars in dwarf 
galaxies by increasing feedback had the trouble of ejecting 
their gas in high redshift star bursts, resulti ng in isolated 
dwarf spheroidal galaxies ([Sawala et al.ll2010l h rather than 
dwarfs with prolonged star formation histories and signifi- 
cant amouiitsof H^gas at z = 0, as observed in the field 
fe.g. lGeha et al.ll2006h . Our feedback implementation is able 
to maintain significant amounts of cold gas in the star form- 
ing region, to regulate star formation whilst driving outflows 
and gas recycling via galactic fountains. 

The match to the baryonic mass-metallicity relation 
(panel 8d) provides added confidence in the baryon cycle 
that our feedback creates. The balance between fresh gas 
accretion, metal injection from star forming regions, metal 
enriched outflows and recycling conspires to result in the 
correct gas metallicities at each mass range. We show in 
Section 5 that this cycle of metals also results in OVI abun- 
dances in the circum-galactic medium (CGM) that match 
observations. Further detailed studies of the metallicities 
and their distributions will be made in forthcoming stud- 
ies, placing greater constraints on our baryon cycle. 

The low feedback case (SG1LF) is also particularly 
interesting. The huge increase in number of stars formed 
means much more mass in the inner region, so the rotation 
velocity (measured at 3.5 scale- lengths) is much higher than 
in the high feedback case. Yet the increased stellar mass also 
results in higher luminosity: the values conspire so that it 
remains on the Tully-Fisher relation. Of course, the peaked 
rotation curve means a there is dependence on the radius 
at which rotation is measured. Beyond the mass-metallicity 
relation, it requires the size-metallicity relation to break the 
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Figure 7. Scaling relations of the simulated galaxies overplotted 
on observational data. All plots are shown on a log scale. Panel 
(a) Stellar mass (M*) against total mas s (M v , r ). The thick red 
line shows the (semi) empirical relation (Sawala et al. 2011 ). All 
other plots use a single observational data set (|Courteau et all 
2007). Red lines are fits to the observational data while dashed 
red lines include 97% of data. Panel (b) Rotational velocity (V c ) 
against Luminosity (Lj) in the I-band (the Tully-Fisher relation). 
Panel (c) Size (S), or disc scale-length, against Lj. Panel (d) S 
against V c . The four fiducial runs are shown as yellow diamonds, 
the low resolution runs as green diamonds and the low feedback 
run as a red triangle. 



degeneracy of the high and low feedback cases both fitting on 
the Tully-Fisher relation. The low feedback run is too con- 
centrated. Thus the shape of the rotation curve and the three 
parameters size-luminiosity-rotation a re required to a ssess 
simul ations. This was pointed out in iDutton &; Courteaul 
(2008). The low feedback run is also deficient in HI relative 
to its R-band luminosity (panel 8c) , due to turning too much 
of its cold gas into stars. 



5 METAL ENRICHED GASEOUS CORONA 

One way to test the extent of outflows, and hence the 
feedback in our model, is to compare the metal enriched 
gas surrounding the s imulated gala xies with observa tions. 
IProchaska et alj (|201lh and iTumlinson et al. | (|201lh find 
that column densities of OVI extend out to 300 kpc from 
star forming galaxies. Their observations extend to galaxies 
as faint as 0.01 L*, making them especially useful to compare 
with our models. Figure [9] shows the surface density profiles 
of OVI gas as a function of the impact parameter, p, from the 
centre of the simulated galaxies. The obser ved OVI column 
densi ties for sub-L* g alaxies at z ~ from ( Prochaska et al.l 
l201lL blue dots) and ([Tumlinson et alj|2011 , green dots) are 
shown overplotted on the model data. Each panel is labeled 
with the V-band luminosity of the simulated galaxies and 
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Figure 8. Simulated galaxies overplotted on observational data. 
Panel (a) Total baryonic mass (Mbaryon ), consisting of all stars 
and c old gas, plotted against V c (observed data from McGaugh 
2005). Panel (b) Sp ecific star formatio n rate against stellar mass 
(observed data from Salim et al. 2007). Panel (c) The ratio of the 
mass of neutral hydrogen gas (H I) to luminosity in the R-band 
(Lr) plotted against the R-band magnitude (observed data from 
IVerheiien fe Sancisi 2001). Red lines are fits to the observational 
data while dashed red lines include 97% of data. Panel (d) The ra- 
tio of Oxygen to Hydrogen . 12+log(Q/H), again st total baryonic 
mass (observed data from Tremo nti et al.ll2004t ). The four fidu- 
cial runs are shown as yellow diamonds, the low resolution runs 
as green diamonds and the low feedback run as a red triangle. 



the virial radius, indicated by the vertical green lines. Our 
feedback implementation produces extended, metal enriched 
gaseous coronae that extend to impact parameters of ~300 
kpc, even around low mass halos, matc hing observed absorp - 
tion line features. We refer the reader to lStinson et al.l ([201 lh 
which presents full details of our comparisons between sim- 
ulations and observations and shows that the low feedback 
case is significantly deficient in OVI within its CGM. Note 
that the affects of metal diffusion on the distribution of met- 
als were explored in IShen et all (|2010h and they find that the 
effect is minimal these results. 



6 CONCLUSIONS 

We use the same physical model, including feedback imple- 
mentation, at the same resolution to simulate four galaxies 
over a wide mass range that match a range of galaxy prop- 
erties and scaling relations. Our model includes radiation 
energy from massive stars, as well as supernova energy. We 
adjusted the star formation efficiency parameter c* so that a 
single simulate d galaxy, SG3 , matched the stellar ma ss-halo 
mass relation ([Moster et all 120101 : iGuo et aljfanoh . With 
the additional energy sources driving outflows from star for- 
mation regions in our updated model, the complex interplay 
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Figure 9. Radial profiles of the column density maps of OVI for 
the 4 simulated galaxies. The large dots are observations of L < 
0.1 galaxies from Prochaska et al. (2011) (blue, top panels) and 
in the lower panels the 0.1 <L<L galaxies from Prochaska et 
al. (2011, blue) and Tumlinson et al. (2011, green) galaxy samples. 
The solid vertical green (red) line represents the virial radius in 
the upper (lower) panels for each of the four halos. 



between cooling gas, star formation, radiation energy injec- 
tion from massive stars, stellar winds and supernova explo- 
sions, outflows, and the recycling of gas (the baryon cycle) 
have been balanced in our model to (i) retain the correct 
amount of baryons, blowing out large amounts of gas from 
the central regions, particularly in low mass galaxies, (ii) 
reproduce star formation rates and histories to attain the 
correct galaxy luminosities, stellar masses and colours for a 
given mass galaxy, (iii) ensure that low mass galaxies are gas 
rich compared to high mass galaxies, relative to their stel- 
lar population, even though significantly more gas is blown 
out of the low mass galaxies, (iv) attain an amount and dis- 
tribution of angular momentum in the star forming gas to 
result in the correct morphology and disc sizes, and (v) cor- 
rectly balance inflows of fresh gas with the recycling of gas 
from star forming regions, in order to correctly match the 
mass-met allicity relation. 

We have previously shown that the out flows in our 
mode l remove low angular momentum gas (iBrook et al. 
l201lh . resolving the "angular mo mentum p roblem" which 
has plagued simulations ( e.g. iNavarro fe Steinmet j [2000; 
IPiontek fe SteTn metz 201 If). Some of the heated gas is 
ejected from the galaxy entirely, particularly at early times 
and from the low mass galaxies, while some of the gas is 
blown only as far as the dark matter halo that surrounds 
the galaxy, and may cool back down to t he star forming re - 
gion at later times and form disc stars ([Brook et al.ll2012h . 
In our four simulated disc galaxies, outflows from the central 
regions also drive an expansion of the central dark matter 



profiles, resulting in flat dark matter central density profiles 
which better match observed galaxies, rather than the very 
steep density profile that are predicted from dark matter 
simulations even in our most massive disc galaxy simula- 
tion. 

By re-calibrating the input star formation efficiency in 
order to attain the same star formation rate in a lower reso- 
lution simulation (SG3LR) as in the fiducial run (SG3) with- 
out changing the feedback implementation, we were able 
to attain a very similar baron cycle and hence a galaxy 
with very similar properties at two resolutions. Our calibra- 
tion for this lower resolution also resulted in an L* galaxy 
(SG5LR) that matches the scaling relations. By contrast, 
a low feedback run (SG1LF) was shown to form too many 
stars, particularly in the central regions. Our study shows 
that it is necessary to examine a number of relations to 
determine the success of simulations, and we recommend 
plotting at minimum the size-luminosity-rotation velocity 
prop erties rather than simply the Tully-Fisher relation (see 
also button fe Court eaull2008l ). 

Our model suggests that large scale outflows of gas from 
the central regions of galaxies are the primary driver of 
galaxy scaling relations, as well as shaping disc morphology. 
The model makes the assumption that radiation energy from 
massive stars regulates star formation and stir up gas in star 
forming regions, allowing SNe remnants to expand and drive 
outflows. We show that the scale of outflows invoked in our 
models matc hes the observed absorption line features of lo - 
cal galaxies ([Prochaska et al ] |201ll : iTumlinson et al.ll20lTh . 
A detailed comparison between the simu lations and these 
observations is made in IStinson et all ([201 lh . where we show 
that lowe r feedback models, mor e akin to those commonly 
used (e.g. lScannapieco et aL|[201lh , are significantly deficient 
in extended OVI, particularly in low mass galaxies. 

Only in a very narrow mass range does the gravita- 
tion of host dark matter halos balance the processes causing 
these outflows, allowing disc galaxies to form. In the lower 
mass halos, the significant energy feedback drives increased 
turbulence, leading to more irregular galaxy morphologies 
(SGI, SG2). 
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